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ONE-DIMENSIONAL NUMERICAL ANALYSIS OF THE TRANSIENT 


RESPONSE OF THERMAL PROTECTION SYSTEMS 

By Robert T. Swann, Claud M. Pittman, 
and James C. Smith 
Langley Research Center 


SUMMARY 


Differential equations governing the transient response of thermal protec- 
tion systems to a hyperthermal environment are presented. These equations are 
expanded into finite-difference equations which are suitable for numerical solu- 
tion. The equations provide for three layers of different materials, the first 
two of which may have moving boundaries. Concentrated heat sinks, such as 
metallic structures, may be located at the back surfaces of the second or third 
layers or of both layers. 

The analysis was developed primarily for charring ablators but is also 
applicable to impregnated ceramic, subliming, and heat-sink thermal protection 
systems. The principal difficulty encountered in numerical analysis of charring 
ablators is the extensive computer time required to obtain solutions. Provision 
is made in the numerical equations to introduce options which reduce computer 
time. The errors resulting from these options under various conditions are 
discussed. Good agreement is obtained between numerical results and exact 
solutions. 


INTRODUCTION 


An adequate thermal protection system may constitute 20 to 50 percent of 
the total reentry weight for vehicles which must enter the earth* s atmosphere 
at super circular velocity. Preliminary studies (ref. l) indicate that charring 
ablators provide the most efficient thermal protection shield for most of the 
vehicle. Other materials may be required for certain vehicle areas. 

Extensive, experimental investigations have been conducted on the perform- 
ance of charring ablators. (See refs. 2 to 6. ) Most of this work is concerned 
with the overall performance of the material in a given environment. The theo- 
retical analysis of charring ablator systems is very complex by comparison with 
the analysis of other heat- shield systems, and many approximations and assump- 
tions must be made to obtain solutions. Several procedures for the numerical 
analysis of the performance of these materials during reentry have been reported. 
These analyses differ primarily in the treatment of thermal decomposition 



processes. In the more rigorous approach (refs. 7 and 8 ), details of the chem- 
ical kinetics of the reactions through the depth of the material are considered. 
In the more simplified approach to the problem (ref. 9)> it is assumed that the 
decomposition processes occur in a single plane at a temperature which may be 
fixed or a function of the rate of decomposition. Since it has been shown that 
the thickness of the decomposition zone is very small compared with the total 
thickness of the material (ref. 10), this simplifying assumption appears to be 
valid. In the present analysis, the simplified approach is followed. 

The most difficult problem encountered in analysis of the performance of 
charring ablators is the formulation of a quantitative expression for the rate 
of char removal. It appears that the char may be removed chemically (oxida- 
tion), thermally (sublimation), mechanically (spalling, aerodynamic shear), or 
by a combination of these methods, depending on the specific material involved. 
Results of an investigation of oxidation as a mechanism in the performance of 
charring ablators are presented in reference 11. A theoretical model which pre- 
dicts spalling of the char, which is observed in some tests, is presented in 
reference 12. It is also ■ anticipated that aerodynamic shear forces have a pro- 
nounced effect on the removal of low-density chars. An analytical program must 
have considerable flexibility if it is to be suitable for investigation of the 
various mechanisms of char removal that may be operative. 

In this paper, equations are derived for calculating the thermal response 
of charring ablators to reentry heating conditions. Sufficient options are pro- 
vided so that not only may all methods of char removal be considered but also, 
by proper manipulation, the response of impregnated ceramics, subliming abla- 
tors, and heat- sink materials to severe heating conditions can be determined. 

The equations have been programed for solution by an electronic data processing 
system, and numerical examples are presented to demonstrate the accuracy of 
various approximations. The equations presented here are similar to those in 
reference 9* The primary difference is that in this paper the finite-difference 
equations are derived for fixed points in a moving coordinate system whereas in 
reference 9 moving points in a fixed coordinate system are used. The present 
paper also provides greater flexibility in the boundary conditions which can be 
imposed. 


SYMBOLS 


The units used for the physical quantities defined in this paper are given 
both in the U.S. Customary Units and in the International System of Units (SI). 
Factors relating the two systems are given in reference 13. 

A constant in mass loss rate equation corresponding to specific 

reaction rate 

B constant in exponential of mass loss rate equation corresponding to 

activation energy 


2 


b 


c e 

Ci+jj 

c w 

Cl, C 2 


thickness of layer at which transfer is made between complete equa- 
tions and boundary equations only (see appendix C) 

oxygen concentration by weight external to boundary layer 

Ci+j+m heat capacity of heat sink (P c p x ) heat sink 
oxygen concentration fay weight at wall 
constants of integration (eq. (5*0) 


C P 

D 

F 

G 

g 

H c 

He 

Hw 

Ah c 

Ahf 

Afap 

i 

A 

K 

k 

k n, n+1 

l 


specific heat 

specific heat of gaseous products of pyrolysis 
parameter defined fay equation (55); ^(^c^p + m P G p) 

heat-generation function 

constant in equation (9d) 

initial temperature distribution 

heat of sublimation of char 

local enthalpy external to boundary layer 

local enthalpy of fluid at wall 

heat of combustion 

heat capacity of coolant 

heat of pyrolysis 

number of stations in char, including one at each surface 
number of stations in uncharred material, including one at faackface 
reaction- rate constant for oxidation of char 
thermal conductivity 

thermal conductivity evaluated at temperature of point midway between 
stations n and n+1 

distance between stations 
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m 


m 

*c 

m p 


number of stations in insulation, including one at back surface 
mass loss rate, a , c m c + a-pip 

rate of char loss 

rate of loss of uncharred material 


m ox 

*Le 

q 

^aero 


rate at vhich oxygen diffuses to surface 

Levis number 

total pressure at vail 

heating rate 

net aerodynamic heating rate at the surface (eq. ( 25 b)) 


convective heating rate to a cold vail vith no mass transfer 
hot-vall convective heating rate corrected for transpiration 
radiant heating rate 
effective nose radius 

unit step function: S = 0, 0 < T; S = 1, 0 ^ T 

temperature 

temperature to vhich back surface radiates 
temperature at finite-difference stations 
temperature of pyrolysis 

T i-fj' T i+j+m temperature at vhich cooling system is activated 


*0 

q G, net 

% 

R 

s(e - t) 

T 

t b 

T n 

Tn 


T 1 

t 

V 


temperature at which char sublimes 
time 

velocity 
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AW* 


y 

z 

a 

a c 

(3 

G 

S 

e 

A 

I 

P 

cr 


rate of coolant consumption 
thickness 

distance from initial outer surface to outer surface of char layer 
distance from initial outer surface of char layer 
distance from back surface of shield 
absorptivity of char surface 

weighting factor for transpiration effectiveness of char mass loss 

weighting factor for transpiration effectiveness of pyrolysis products 

either 0 or 1 depending on whether transpiration or ablation theory is 
used 

emissivity 

transformed coordinate for uncharred material (eq. ( 23 b)) 

transpiration coefficient 

temperature in differential equations 

weight of char removed per unit weight of oxygen 

transformed coordinate for char layer (eq. ( 23 a)) 

density 

Stefan-Boltzmann constant 


Subscripts: 

o initial value 

i, j,m number of stations in first, second, and third layers (see fig. 2) 
n integer 

s stagnation point 

oo free stream 

CALC calculated by numerical method 

EXACT exact solution 
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Unprimed symbols refer to char layer unless otherwise 
primes refer to uncharred material; double primes refer to 


specified; single 
insulation. 


ANALYSIS 


The thermal protection system that is to be analyzed is shown schematically 
in figure 1. Although this discussion is confined to a charring ablator system, 

all the concepts and equations apply equally as well to any 
other thermal protection system composed of not more than 
three primary layers . For a charring ablator system, the 
outer (heated) layer is the char, the center layer contains 
the uncharred material, and the third layer consists of 
insulation. Heat sinks can be located at the back of the 
second or third layers or at both locations. 

Char Layer 


_ WHorChar_ Surface^ . ^ 


Uncharred Material 


The outer (char) surface is subjected to aerodynamic 
heating. The char layer provides both insulation and a 
high-temperature outer surface for reradiation. The heat 
passing through this layer is partially absorbed by pyrol- 
ysis at the interface between the char layer and the 
un charred material, and the remaining heat is conducted into 
the uncharred material. The gases generated by pyrolysis 
transpire through the char layer and are injected into the 
boundary layer. The gases are heated as they pass through 
the char, and this heat removal from the char layer reduces 
the quantity of heat conducted to the pyrolysis interface. 
When these gases are injected into the boundary layer, the 
convective heat transfer is reduced. This reduction in con- 
vective heating is the same effect as that obtained with 
simple subliming ablators. In addition to the gases produced by pyrolysis, the 
carbonaceous residue remaining at the interface adds to the thickness of the 
char layer. While the processes of pyrolysis, transpiration, and injection are 
underway, char removal may also be taking place as a result of thermal, chemical, 
or mechanical processes. Thus the total char thickness may increase or decrease 
depending on the relative rates of formation and removal of the char. The var- 
ious processes discussed are related quantitatively in the following sections. 


Figure 1.- Schematic 
diagram of system 
employing charring 
ablator. 


Differential Equations 


It is assumed that thermal properties in a given layer of material are 
functions only of temperature, that all heat flow is normal to the surface, and 
that gases transpiring through the char are at the same temperature as the char. 
Then the governing differential equations (from ref. 9 ) for the char layer 
(x = y ^ x + x) are as follows: 
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dy\ 


- Se 
Sy 


m p c p 


pep 


Se 

5 t 


Heat conducted 


Heat absorbed by 
transpiring gases 


Heat generated Heat stored 


for the uncharred layer 


x + x 


y ^ * c 


+ x i)^ 


-A/k< Ml.) 
Sy / 


p 'c 


, Ml 

p at 


and for a layer of 


insulation 


(*o + x o = ^ = x o + x o + x o)> 


_a 


M 



P"< 


Mil 

St 


( 1 ) 


( 2 ) 


(3) 


The thicknesses of the layers to which the first two of these equations apply 
vary with time in a manner which is determined by the boundary conditions. 


Initial Conditions 

The initial temperature distribution is assumed to be given as a function 
of position: 


e(y,o) = g(y) (4) 

The initial mass-transfer rates must also be specified. It should be noted that 
these values can be other than zero for some cases. 


Surface Boundary Conditions 

Two conditions must be specified at the heated surface. One must determine 
either the rate of removal of material at the surface or the temperature of the 
surface; the other is provided by the energy balance. 

Surface ablation .- In general, the relative importance of the mechanisms 
involved in char removal from specific materials is not well established at this 
time. It has been established, however, that oxidation of the char surface is 
one important mechanism. Spalling of the char as a result of internal pressure 
is observed in some cases. Ablation at a given temperature (that is, sublima- 
tion) occurs if the heating rate is sufficiently high. Ablation of the surface 
may also occur as a result of aerodynamic shear stresses. 

To provide maximum flexibility, provision is made for the following mecha- 
nisms of surface erosion: 

(l) Ablation at a given temperature which may be a function of ablation 
rate (sublimation) 
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(2) Removal of char at a rate which is a given function of time (spalling, 

aerodynamic shear) 

(3) Removal of char at such a rate that the char thickness is a given func- 

tion of time (spalling, aerodynamic shear) 

(b) Ablation as a result of a chemical process (oxidation) 

For ablation at a given temperature, two cases are considered. In one case, 
ablation occurs at a fixed temperature. In the other case, the char mass loss 
rate is an exponential function of the surface temperature. For ablation at a 
fixed temperature, no surface erosion occurs if the calculated surface tempera- 
ture is less than the specified ablation temperature. If the calculated surface 
temperature is higher than the ablation temperature, ablation occurs at a rate 
sufficient to reduce the temperature to the ablation temperature; that is, 
is equal to zero for Tq < Tq, and me is calculated from an energy balance at 
the surface for Tq = Tq. In the second case, char mass loss rate and surface 
temperature are related as follows: 

m c = Ae~ B / Tl (5) 

An equation of this form has some physical significance, because decomposition 
reactions proceed more rapidly at higher temperatures. By an appropriate selec- 
tion of A and B, equation ( 5 ) yields results similar to those obtained by 
specifying an ablation temperature. 

If the rate of char removal is given function of time 

m c = f(t) (6) 

then the rate of char removal is obtained from the input data and the surface 
temperature is calculated from an energy balance. Such a relation might be used 
to compare calculated and experimental results when a more basic quantitative 
relation for the experimental rate of char removal is not available. 

If the char thickness is a given function of time 

x = f(t) (7) 

then the rate of char removal is calculated from this relation together with the 
rate of char formation which is calculated from the conditions at the pyrolysis 
interface. This condition can be used when it is desirable to perform calcula- 
tions for applications in which the char thickness is known as a function of 
time even though mechanisms of char removal may be present which cannot be 
expressed quantitatively. 

It has been shown experimentally that oxidation is an important mechanism 
of char removal. (See ref. 11.) For a half-order reaction, the rate of oxida- 
tion of carbon can be determined from the following equation (ref. lb): 
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m c = Ae" B / Tl ^C w p w 


(8) 


The pressure at the wall must he specified for subsonic and supersonic flow. 
However, in hypersonic flow, the pressure can be related to the stagnation 
heating rate and enthalpy. The stagnation pressure in hypersonic flow is 
approximately (see ref. 15): 


Further (from ref. l6), 


and 


r, _ n n 

£*w, s 12 


V 2 
00 00 



h e 


OC 


Then, P v ^ s can be expressed by the relation 



( 9 a) 

( 913 ) 


( 9 c) 


( 9 d) 


where the constant of proportionality is 


G = 410.72 


ft?. 


sec 


lb £ 


-atm _ ^ 2 m^-sec 2 -atm 


kg^ 


The pressure at the wall is therefore given by 


P w 



( 9 e) 


where depends on vehicle attitude and body location. The rate at which 

■^WjS 

char is removed by oxygen must be proportional to the net rate at which oxygen 
diffuses to the surface. From reference 17 this rate is 


%x - 


w 0.6 

^Le q 0,net / G 
he “ hw e 



( 10a) 


As shown in a subsequent section, q^ ne ^. 
rate corrected for transpiration (see eq. 


is the hot-wall convective heating 
(13)); that is, 
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(10b) 



By eliminating the concentration of oxygen at the wall C w in equa- 
tions (8) and (10a), the rate of removal of char by oxidation is found to be 


m„ 





1 (h e - h w )K 2 p w ^ 

Rh e - h w )K 2 p w 

2 

+ ^ 2 p w c e 

f 0.6 1 

^0-6 

[ Oc.net^Le f 

q C,net^Le 

J 


(10c) 


where K = Ae 




The equation for a first-order oxidation reaction is obtained similarly. 
The resulting equation is 


K Pw C e 

Kp w (h e - h w ) 

^ 0.6 

Pc,net^Le 


(lOd) 


Surface location .- "When char removal occurs, the char surface moves with 
respect to a coordinate system fixed in the material. The distance between the 
surface of the char and the initial surface location is given by 


x = 



(id 


The thickness of the char at any time is equal to the initial char thickness, 
plus the thickness of char formed by pyrolysis, less the thickness of char 
removed; that is. 


x = x Q + 




( 12 ) 


Surf ace energy balance . - The heat input consists of convective and radiant 
heating. This energy must be accommodated at the surface by a combination of 
f our mechanisms : 
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(1) Blocking by mass transfer into the boundary layer 

(2) Reradiation or reflection from the surface 

(3) Conduction into the material 

( 4 ) Sublimation of the char 

The effect of mass transfer on heat transfer has been studied extensively. 
With low-mas s-transfer rates it is found that the reduction in heat-transfer 
rate is directly proportional to the product of the mas s-transfer rate and the 
enthalpy difference across the boundary layer. With high-mass-transfer rates, 
which may occur when a large fraction of the heat input is radiant, the linear 
approximation is no longer adequate and it is necessary to use a higher order 
approximation. A second-degree approximation is derived in appendix A. 

The surface energy balance, expressed in a form in which either approxima- 
tion to the blocking effectiveness can be selected, is as follows: 




Cold vail 
convective 
heating rate 



- (1 - 3) 


0.724 + <V*p) - °- 15 ^ ( a o A c + V’p) 2 ' ?l(v= + 


Aerodynamic blocking 


Net convective heating rate 




+ 


aq R 


+ [1 - s(e - t-J] Ah c = 

v ^ ' 


L 

a€ 1®1 


Radiative Combust ive Reradiation 

heating rate heating rate 



s(e - t^AcHc 
' v ' 


Conduction Heat of sublimation 

to interior of char 


(13) 


If transpiration theory (second-degree approximation, appendix A) is used, 

P is equal to zero. For linear ablation theory, P is equal to 1. In either 
case, the heat absorbed by vaporization of the char H c and the heat of com- 
bustion of the char Ah c are considered separately. The coefficients a, c and 
ap can be used to differentiate between the blocking effectiveness of the gases 
produced at the surface and at the pyrolysis interface. Evaluation of these 
coefficients is discussed briefly in appendix A. 

The heat transfer to the outer surface is assumed to be a given function 
of time and consists of the cold-wall convective heating rate q^ and the radi- 
ant heating rate q R incident on the surface. These two components must be 

specified separately because mass transfer at the surface blocks part of the 
aerodynamic heating but, in general, has no effect on radiant heating. Addi- 
tional terms can easily be included in equation ( 13 ) to account for other 
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phenomena which may affect the heat balance at the char surface. For example, 
reference 18 discusses a gas-phase combustion in the boundary layer involving 
the gases of pyrolysis. This effect has not been clearly identified at the 
Langley Research Center and is, therefore, not included in the equation. How- 
ever, phenomena such as this may be important in some cases and their existence 
should certainly be considered. 

Equation (lj) is normally used in this analysis as the boundary condition 
on the temperature at the outer surface. However, when 9 is equal to the sub- 
limation temperature Tq, the specified sublimation temperature provides the 
boundary condition on the temperature and equation ( 13 ) is used to calculate the 
rate of ablation m c . 


Pyrolysis-Interface Boundary Condition 

Energy balance .- The heat conducted to the pyrolysis interface must be 
either absorbed by pyrolysis reactions or conducted into the uncharred material; 
that is, at y = x + x. 


■k^ = 

oy 


m-n Atir 


be' 

dy 


(14) 


In addition, the temperatures in the char and in the uncharred material must be 
equal at the interface; that is, at y = x + x, 

e = 0’ (15) 

Pyrolysis rate . - Two approaches are available for calculating the rate of 
pyrolysis. In the first approach, it is assumed that pyrolysis occurs at a 
given temperature T-j_. If 


®x+x 


< T i 


(l6a) 


then 


m p = 0 


( 16 b) 


If 


®5c+x - T i 

the temperature is known and equation (l4)- is used to calculate the rate of 
pyrolysis. 

In an alternate approach, it is assumed that the rate of pyrolysis is a 
known function of temperature, for example 
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mp = A 1 e 


- B V e x+x 


(IT) 


when 


9 x+x < T i 

In this case, equations (l4) and (17) are solved for both temperature and pyrol- 
ysis rate. The value of T^ is still specif ied, and if this temperature is 
reached, the pyrolysis rate is determined only from equation (l4) so that this 
temperature is not exceeded. 

Pyrolysis-interface location .- As pyrolysis occurs, the interface "between 
the char layer and the uncharred material moves with respect to a fixed coor- 
dinate. Its distance from the initial char surface location is 


x + x = x 0 + / — dt 

Jo p - p 

The instantaneous thickness of the uncharred material is 



(18) 


(19) 


Boundary Conditions at Back Surface of Ablation Material 

A number of conditions can be imposed at the back surface of the ablation 
material, depending on whether additional insulation is provided or some pro- 
vision is made for temperature control. Whether insulation is used or not, the 
ablation material may be attached to a thermally thin plate which functions as 
a concentrated heat sink. 


Three-layer system .- If insulation is used, the temperature of the ablation 
material is equal to the temperature of the insulation at their interface. 


9x o +x i 0x o +x i 

From an energy balance at the back surface of the ablation material. 


-k 


, del „ 
Sy 


= c 


ae' 
i+ J at 


... ae_ 

dy 


(20a) 


(20b) 


Two-layer system .- If no insulating layer is used, the back surface can 
be assumed to be perfectly insulated, cooled at a given temperature, or may 
exchange radiation with a sink of known temperature in the interior of the 
structure. An energy balance yields the following equation: 


I 
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-k' 


de* _ P 

Sy ~ ° i+ j dy 


+ s (e 1 - 



( 21 ) 


The temperature at which the cooling system is activated is T^ + j. The choice 

of conditions is accomplished by making the inapplicable terms equal to zero 
(that is, C 1+ j = 0 and/or T i+ j > 0 f and/or e i+ j = Oj. 


Boundary Condition at Back Surface 
of Insulating Material 

If an insulating material is used behind the ablating material, the 
boundary condition at the back surface (y = x Q + + x”) is similar to 

equation (2l); that is, 

-k" = C i+J+m + s (e" - T i+ j +m )w f Ahf + ae i+J+m [(0") 4 - T 4 ] (22) 


Transformation of Coordinates 

The equations derived in the preceding discussion are similar to those pre- 
sented in reference 9- In reference 9, these equations are expressed in finite- 
difference form and solved in a fixed coordinate system. To maintain a fixed 
number of stations in layers of varying thicknesses, it is necessary to change 
the locations of the stations and to interpolate to determine the temperatures 
at the new locations after each step in the calculation. This procedure not 
only increases the time required to perform the computations, but also intro- 
duces a small error In each step of the calculation. This difficulty can be 
eliminated by transforming the equations to a coordinate system in which the 
finite-difference stations remain fixed, and the coordinates themselves move 
to accommodate the changes in the locations of the surfaces of the different 
materials. 

The y-coordinate can be transformed to £- and ^-coordinates in the char 
and uncharred layers, respectively, by using the following equations: 



( 23 a) 


( 23 b) 


lb 



In this coordinate system the outer surface remains fixed at £ = 0. The inter- 
face is located at £ = 1 in the char and at £ = 0 in the uncharred material. 
The hack surface of the uncharred material is located at £ = 1. A number of 
advantages result from the use of this double transformation. First, the char 
always extends from £ = 0 to £ = 1. Therefore, the temperatures tend to be 
more nearly steady state than would be the case with a coordinate system fixed 
at the surface only. A second advantage is the positive location of the pyrol- 
ysis interface. A similar transformation would also be very beneficial in 
locating the center of the reaction zone when the pyrolysis reactions are con- 
sidered in detail. Because the reaction zone is typically very thin, a very 
fine finite-difference network is required to analyze it. With transformations 
similar to those here, the center of the reaction zone can be located, and the 
fine network can be restricted to this region rather than covering the entire 
range of possible reaction-zone locations. 

In the transformed coordinate system, equations (l) and (2) are as follows 
(for the char layer and uncharred layer, respectively): 


4^-fk & 


x2 5| \ 5f 


(x'Y 


, 1 50 
x hi 


m c c p 


m P c p 


+ I 


m pP 
P 1 - P 


me) Cp 


50 

St 




F = P<=p 


58' 

St 


(24a) 

(24b) 


The boundary conditions are as follows: 


At |=0, 


= ae l 04 “ x + A c 


aero 


s(e - T x ) (h c + Ah c ) - Ah c 


where 


‘Wo = ac lR + “ (1 “ P)|o-T24 ^-(a-c & c + “p^p) 

- 0 . (acme + a pmp)J “ P 1 ! (ac&c + a p%>)|^j 


at i = 1, £ = 0, 

and 


0=0' 


_k 2a. = Ah-n 
x P P 


h 1 58 
x’ 5£ 


( 25 a) 


( 25 b) 

( 26 a) 

( 26 b) 
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and at £ = 1 for only two layers, the condition at the hack surface 
(y = x 0 + x^) is 


k 1 8e' _ aei 

T ~ ^ 


+ S 


( e ' - vj )^ + « i +J [( 9 ’) 4 - 4 ] 


( 27 ) 


x' % i+ j at 

When three layers are used, the conditions at the hack of the second layer are 

0 f =0" (28a) 

and 


ki Mi = n Mi v« Ml 

“x* ^ at ' 8y 


( 28 b) 


If three layers are used, the condition at the hack surface (y = x Q + Xq + Xq) 


is 


-k" = C n - 


80 


8y 1+ J +m at 


+ s^e" - T i+ j +m ^ f Ah f + - TgJ (29) 


Finite-Difference Equations 

The locations of the finite-difference stations are shown in figure 2. 
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Figure 2.- Location of finite- 
difference stations. 


The methods used to change the differential equations to finite-difference form 
are given in appendix B. A summary of the finite-difference equations obtained 
by these methods is presented here. 

The equation at the first station (eq. (B25)) is 


T 2 - Ti 


% 2 T 


+ ^ero ~ ae 


1*1 - *c[s(Tl - T x )(h c + Ah c ) - Ah c ] 


+ ^(" 11T 1 + 18t 2 - 9T 3 + 2%)(m c c p + mpCp) + \ F = pc p \ ^ (30) 

The equations for the first layer (eqs. (B9) and (BIO)) are 


n-l, n Tn ~y~ T - - n+1 ~ + f ( 6T n+l “ 3 T n - 2?n-l - T n+2 ) 


T n " T n+1 . l. 


m c c p + nipCp 


n - If m p p . \ 

+ “ mc J c p 


and 


+ ZF = Zpc 


AT 


n 


P At 


(2 ts n ^ i - 2) ( 31 ) 


T n-1 " T n 


LX — -L- XX -i 

c n-l, n ^ K n, n+1 


T n “ ^n+1 , 1 


+ J(2Tn+l + 3T n - 6T n .q + T n . 2 ) 


m c Cp + mp Cp 


,n-l V . ) 

+ —1^7 “ -c e p 


AT. 


n 


+ ZF - PC P Z At 
The equation at the pyrolysis interface (eq. (B30)) is 


(n = 1 - 1) (32) 


T. ,, - T, 

P( lr ' P V 

K i, i+1 z < *1-1,1 


T i - * 1-1 


+ IF + c-n + 


C T)P 


“P “p + pTTTp 


i(llTi - l8T i _ 1 


+ 9T ± _ 2 - 2T 1 _ 3 )J + ^[K-HTi + 18T 1+1 - 9T i+2 + 2T 1+5 ) 

- 2Ahp^ = ^pCpZ + p'cpZ'j 


vATt 

/At 


(33) 
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The equations for the second layer (eqs. (B15) and (Bl 6 )) are 


k’ 


■^n -1 “ ^n , t ^n ~ -*-n+l 


- k* 


n- 1 , n 1 1 n, n +1 2 1 


• . t 


and 


' ^n -1 “ ^n ^ 1 


n-l,n 2 . 

§(’ 


n, n +1 i . 


+ ^(2T n+1 + 


1 J 

p' - P 


(i + 1 ^ 

T n +1 


i + j - 

n V’ C P 


= 2 'p’Cp 


AT r 

AT 


•zm 6T -, + T X 1 + J ~ n V’S _ » ^n 

3T n - bTn.! + T n _ 2 j . p' _ p 1 p c p £t~ 

(n = i + j - 1 ) ( 35 ) 


The equations at the hack of the second layer are for two layers 
(eq. (B33)) 

k' T i+j-l ~ T i+j Lk JA , 

k i+j-l,i+j z x ~ 0 e i+j \™i+j “ T B/ " S ( T i+j 


( T i+j - T i +J )«f 

S p c p T 


+ c 


AT. . . 
J-+J 

At 


and, for three layers (eq. (B 36 )) 


HP . rn 

k ' i+j-1 X i+J v " 

^i+j- 1 , i+j ” K i+j,.i+j+l 


T i+j ~ T i+j+l 


, it 1' ^ t. it l" n \^i+j 

-lP c p 2 + P c p 2 + c i+j) At 


(36) 


(37) 


The equation for the third layer (eq. (B19)) 


, •* -^n-1 ^n ^n “ ^n+1 _ n m, » ^®n 

n-l,n p k n, n +1 " p V At 


(38) 
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The equation at the back surface of the third layer (eq. (B39)) is 


i+j+m-1, i+j+m 


Tj+j+m-1 ~ T i+j+m 


cr£ i+j+m (^i+j+m ~ t b) + s( T i+j+m “ T i+j+m)^ iW f ^f + C P ~ + C i+0+mj 


N^i+j+m 

y~At 

(39) 


In some cases it is desirable to use the so-called thin-layer options, 
presented in appendix C, because of the considerable saving of computer time 
involved in their use. When these options are used, the layer in question is 
divided into two half-elements and the interior stations in the layer are no 
longer considered. The first layer is considered "thin" when x is less than 
some prescribed thickness b, and the second layer is considered "thin" when 
x* is less than some prescribed thickness b ! . The finite-difference equa- 
tions used with the thin- layer options are as follows. 

At the front surface, method 1 gives the following equation when x ^ b 
(eq. (Cl)) 


T i - T 1 4 

k l,i x + ^aero " CT£ 1 T 1 “ m c 


! ( T 1 " T l)( 


H c + £h c ) - Ah c 


T i " T 1 /. . . N x x AT 1 

+ J (" C °P + “P^) + 2 F = PC P 2 At" 


(40) 


The equations at the interface are, for x ^ b and x' > b' (eq. (C2)) 


*1,1+1 ^ 1 T- T -- ~ 2kl,i + ^ + “p{( £ P + p^-pJ(T i - T l) 


cAp’ 


P* - P 


|(- 11T 1 + 18t 1+1 - 9T i+ 2 + £T 1+J jl - 22h p | - (pc^ + p'o^l')^ 


(41) 


for x -g b and x'^b’ (eq. (C3)) 


2k. 


T i+j " T i 


Ti - Tq 


“ n 

1^ 1 X 


+ xF + mp 


P C T 


1, i+ j x 

- T i) - a^pj - (p=p* + p V) 


C P ^ ~ 


P “ P, 


(Ti - Tq) 


ATq 
At 


(42) 
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and for x > b and x* < b' (eq. (C4)) 


2k 


i + j 


T i+«j " T i 


- 2k, 


- Ti_! 


x 


+ IF 


||(llTl - 18T 1 _ 1 + 9TJ..2 - ZT^j) 

p ^i-(T 1+3 - Ti) - aap'j) = (pc p i + p'c^'jlr 


P C P 

+ Mi d p 


(^3) 


The equation at the hack surface of the second layer for two layers 
(eq. (C5) ) is 


* — / k k \ / — 

k i,i+j — p 1 - “ e l + j) T i+j - t b] " s ( T i+j - T i+j)»f 


m ('pV c . V'h+J 

\ p C P 2 Ul+J ) At 


(44) 


and for three layers (eq. (C6)) is 


%i+j 


T i ~ T i+j , " 


i+j^ i + j + l 


T i+j ~ Tj+j-i 


P ,e i T + ^ T + c i + 3 ) 




At 


(45) 


Calculation of Temperatures at Fixed Points 

It is possible to calculate the temperature at any number of fixed points 
Z n where Z is the distance from the back surface. (See fig. 2.) The tem- 
perature at is calculated as follows (by using linear interpolation): 

When 


0 < Z n ^ x" 


find station N such that 


Then 


(i + j + m - N)z" ^Zn^Ci + j+ m- N- l)z" 


T( z n) = t N + (% " %+l) 


z n " (•*- + J + m “ N)z" 


(46) 
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When 


x" ^ Z n < x* + x" 

find N such that 

x" + (i + J - N)z' ^ Zn > x" + (i + o - N - 1)2’ 

Then 

Z n - x" - (i + J - N)l 1 


T ( z n) = % + (T N - T n+1 )- 


l 


(47) 


When 


x T + x" ^ Z n ^ x + x l + x" 


find N such that 

x' + x" + (i - N)z ^ Z n ^ x' + x" + (1 - N - 1)Z 

Then 

/ \ / \Z n - X* - x" - (i - w)z .. a . 

T(z n ) = T N + (T N - T n+1 )— - (48) 

If the outside surface moves past the largest Z n because of ablation of the 
surface, then the temperature at this point is meaningless and is, therefore, 
not calculated. 


RESULTS AMD DISCUSSION 


Comparisons have been made of the numerical results obtained when the 
various options presented in this report are used. Comparisons with exact solu- 
tions have also been made. The results of these studies are presented and dis- 
cussed in the following sections. 


Comparison With Exact Solutions 

The accuracy of the numerical analysis is determined for two cases for 
which exact solutions are available. One case provides an exact solution for 
a homogeneous nonablating heat-sink material subjected to a suddenly applied 
constant heating rate. The other case is for quasi- steady- state ablation. In 
both cases, constant material property values and environmental values are used 
and re radiation from the surface is not considered. 
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Heat-sink case . - The exact solution for the heat-sink case is obtained 
from equation (All) of reference 19. When written with the symbols used in 
this paper, the equation is 


T 


(if? 



- To = 


q.(x + x' ) 


kt 


|pc p (x + x')‘ 



X + X 1 / 


1 

6 



(-l) n 


cos 


rut 1 


X + X 


exp 


2 2 
n^rt^ 


kt 


pCp(x + X 


1 


(49) 


The values of q and material properties used to obtain a solution to equa- 
tion (4-9) are given in table I. These values do not represent a real material 
but were chosen to simplify the calculation. 


TABLE I.- VALUES USED FOR HEAT-SINK EXACT SOLUTION 


Quantity 

Value in 
U.S. customary 
units 

<1 

100 Btu/ft 2 -sec 

t . . 

0.1 ft 

k 

0.01 Btu/ft-sec-°R 

P 

10 lb/ft3 

C P 

0.1 Btu/lb-°R 

To 

530° R 


Value in 
SI units 


1135 kw/m 2 
3 . 048 cm 
62. 4 w/m-°K 

l6o kg/m^ 
0.4l84 kj/kg-°K 

294° K 


When a comparison is made between surface temperature histories calculated 
from this equation and calculations made with the finite-difference equations, 

y 

the following results are obtained. At the heated surface = 0, the 

x + x 

error is about -0.4° R (-0.2° K) out of about 1863° R (1035° K) when T 0 is 
530° R (294° K). 

Figure 3 shows a comparison between the finite-difference solutions and 
the exact solution for the first few time steps of the calculations. The cal- 
culations shown in figure 3( a ) use a time interval of l/k09 6 second. It can 
be seen that after about 1 6 time steps (l/25 6 second each), the agreement is 
adequate . 

The results obtained when the time interval is cut in half are shown in 
figure 3(h). It can be seen that the error in surface temperature is reduced 
considerably for the same heating period. 
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(a) At = 1/^096 second. 



(b) At = 1/8192 second. 

Figure 3*- Comparison of surface temperatures for 
initial time steps. 


Quasi- s teady-state ablation case ,- A quasi-steady-state case is one in 
which the front surface and pyrolysis interface recede at the same rate} thus, 
the char thickness is constant. If conditions are also chosen such that there 
is no heat transfer into the uncharred layer, an exact solution can be obtained. 


The transformed heat balance equation for the first layer is (eq. (24a)) 

1 5 


L Se\ , 1 Se 

A .. A ,/ ¥ .\ 

c Sr] * sr 

m C Cp -r IUpCp "r T p " m c/ c p 


, 1-1 W w 

+ F = pep — 


Se 

St 
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With quasi-steady-state ablation, x, nip, and m c are constant. 

- m c = o] 


Further, 


n>pP 


P' - P 


„ o 


St 


J 


Therefore equation (2ka) becomes (with F = 0): 

i a L ae\ , 1 ae/. „ , • - \ _ n 

^2 dFr dI7 5c iFv p p ) " 

Then with constant k, Cp and Cp, equation ( 51 ) can be written 


afe + D ffi.o 
at 2 at 


where 


D = l(“c c p + Vp) 

The solution of equation ( 52 ) is (0 being replaced by T) 

T = Cp + C 2 e~ Z ^ 

and Cp and C 2 can be found from the boundary conditions: 

T = T c 


T = T-r 


(i = 0) 
(I = 1) 


(50) 


(5D 


(52) 


(53) 


(5^) 


The final solution is 


T = Tp - T c e" P + (Tq - T p )e 
1 - e- p 


-D| 


(55) 


To determine the mass loss rates m c and Ap, and the char thickness x 
associated with this equation, the following equations are used: 
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p 


m p p 

P ! - 


= m r 




tfje-o = ^( Mc " q —)} 

del = _ ^p x 

S|J| = i k 


(56) 


In the last equation, it is assumed that — 


Z=0 


is negligibly small; that is. 


there is no heat flow into the second layer. By solving equations ( 55 ) an & 
( 5 6 ) for m c , Ap, and x, the following equations are obtained: 


m c 



^•aero 




+ 



Ap 





PH C 

P* - P 


'N 


> (57) 



The values of laero anc ^ material properties used to obtain a solution 
to equation (55) (by using eqs. (53) an <3. (57)) ara given in table II. These 
values do not represent a real material but were chosen to simplify the 
calculation. 


TABLE II.- VALUES USED FOR QUASI - STEADY - STATE EXACT SOLUTION 


Quantity 

Value In 

U . S . Customary Unit s 

Value in 
SI Units 

q. 

69*3 Btu/ft^-sec 

785 kW/m 2 

k 

1.0 X 10 rh Btu/ft-sec-°R 

0.624 w/m-°K 

P 

20 Ib/ft^ 

320 kg/m^ 

P* 

40 lb/ft^ 

64 0 kg/m^ 

C P 

0.5 Btu/lb-°R 

2.09 kj/kg-°K 

fp 

0.5 Btu/lb -°R 

2.09 kj/kg-°K 

Tl 

4000° R 

2222 °K 

T i 

1000° R 

556 °K 

Ah p 

1000 Btu/lb 

2.324 Mj/kg 

He 

1000 Btu/lb 

2.324 Mj/kg 


25 




X 


(a) Mass loss rate. 


The results of a comparison 
between the quasi -steady- state exact 
solution and the numerical analysis 
are shown in figure 4. For the case 

p 

considered, the ratio — - = 1 

p 1 - p 

so that m c = nip. Figure 4(a) shows 

the ratio of the mass loss rate 
obtained from the numerical solution 
to the value obtained from the exact 
solution plotted as a function of the 
distance l between stations. In 
the first material, the error is less 
than 0.7 percent even with large 2. 



QQfiS ' * •— 1 1 J 

• yyO:> 0 .05 .10 .15 .20 .25 

x 

( b ) Temperature . 

Figure 4.- Error in quasi- steady- state solu- 
tion (x = 0.01 ft = 0.3048 cm). 


The ratio of the temperatures 
obtained from the numerical solution 
to those obtained from the exact solu- 
tion plotted as a function of the dis- 
tance between stations in the first 
material is shown in figure 4(b) for 
three locations through the thickness. 
The error is greatest at the center 
location, as would be expected, since 
the temperature is fixed at the sur- 
face and interface. Note that the 
error is less than 0.2 percent even 
for large 2 . 


Thin- Layer Options 

The time interval for which the numerical solution is stable is propor- 
tional to the square of the distance between finite-difference stations. (See 
appendix C.) In general, with the present formulation for the first deriva- 
tives, the distance between stations is less than or equal to one-fourth the 
thickness of the layers. In many cases, the initial value of x, or the final 
value of x 1 , is very small and excessive computer time is required. To reduce 
computer time in such cases, equations which eliminate the interior stations 
when the thickness of a layer is less than a specified value are provided. 

These equations are derived in appendix C. Specifically, when x ^ b (referred 
to as the b-option), the interior stations in the first layer are not considered 
and the entire layer becomes two half elements of thickness x/2. If the first 
layer reaches a reasonable thickness due to pyrolysis at the interface, the 
thin-layer option is no longer appropriate and the program reverts to the stand- 
ard equations when this value of b is exceeded. A b-value, which is one 
order of magnitude larger than the initial value of x, generally provides 
acceptable accuracy and significantly reduces computer time. Similarly, when 
x ! ^ b 1 (b '-option), the interior stations of the second layer are discarded. 

A b'-value, one order of magnitude less than the initial value of x ! , is gen- 
erally acceptable. 
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Results obtained with the thin-layer options are compared with the standard 
solution in figures 5 and 6 for a typical charring ablator subjected to a con- 
stant heating rate. The values used in the calculations are given in table III. 
These values are representative of a low-density phenolic-nylon charring 
ablator. 



(a) Char layer. 



(b) Uneharred layer. 

Figure 5*- Error in calculated thickness vhen thin- char- 
layer option is used. 
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TABLE III— VALUES USED FOR THIN- LAYER OPTION CALCULATIONS 


Quantity 


1c 

k for - 

500° R (27 8° K). 
1000° R (556° K) 
1500° R (834° K) 
2000° R (1110° K) 
7000° R (389O 0 K) 
C P 





P • • 
P* • • 

Ahp . 

A . . 

B . . 
C e • • 

A . . 

P w , atm 

A 1 . . 
B* . . , 


Value in 

U.S. Customary Units 

100 Btu/ft 2 -sec 

5 X 10~ 6 Btu/ft 2 -sec-°R 
2 x 10"5 Btu/ft 2 -sec-°R 
5 x 10" 5 Btu/ft 2 -sec-°R 
1.0 X 10-^ Btu/ft 2 -sec-°R 
1.0 X 10' 5 Btu/ft 2 -sec-°R 
0.5 Btu/lb-°R 

0.5 Btu/lb-°R 
0.32 Btu/lb-°R 

3.06 X 10"5 Btu/ft-sec-°R 
15 lb/ft 3 

36 lb/ft 5 

250 Btu/lb 
6. 73 x 10 8 lb/ft^sec-atm 1 / 2 

3.9877 x 10^ °R 
0.23 
0.75 
1-5 

1.585 X 10 8 lb /ft 2 - sec- atm 1 / 2 
2.321 X 10^ °R 


Value in St Units 

1135 kW/m 2 

0.0312 w/m-°K 
0.1248 W/m-°K 
0.312 w/m-°K 
0.624 w/m-°K 
6.24 w/m-°K 
2.09 kj/kg-°K 

2.09 kj/kg-°K 
1.34 kj/kg-°K 
0.191 w/a-°K 

24 0 kg/m 5 
576 kg/m 5 
O.58I Mj/kg 

3.29 x 109 kg/m 2 -sec-at m 1 / 2 
2.216 X 10^ °K 


7.74 x 10 6 kg/m 2 -sec-atm 1 / 2 
1.29 x 10^ °K 


i.o6r 


Figure 5(a) shows the ratio of 
char- layer thickness using the thin- 
char- layer option to the char-layer 
thickness using the standard equations. 

As can he seen from the figure, when 
x < b throughout the calculation, there 
is an appreciable error. However, for 
both calculations, the char thickness is 
monotonically increasing and if a b-value 
is used such that x becomes greater 
than b during the calculation, the char 
thickness immediately begins to approach 
the thickness obtained when the option is 
not used. Thus, since most of the com- 
puter time necessary for a calculation 
is usually consumed in the beginning of 
a calculation when the char thickness is 
small, the use of an intermediate b-value 
usually produces a considerable saving of 
computer time without any appreciable 
error. 



(a) Char layer. 



Figure 5(b) shows the ratio of the 
thickness of the un charred layer using 
the thin- charred- layer option to the 
uncharred layer thickness using the 
standard equations. The error which 
occurs when the b-option is used is sig- 
nificant. However, as with the char layer, 
reduces this error considerably. 


(b) Uncharred layer. 

Figure 6.- Error in calculated thick- 
ness when thin-un char red- layer 
option is used. 


the use of an intermediate b-value 


Figure 6 shows the results obtained by using the b ! -option. The value of 
b f used, 0.02 inch ( 0.05 cm), is such that x T becomes less than b 1 at about 
130 seconds. The error shown in the figure can be reduced, if necessary, by 
using a smaller b 1 -value. 


CONCLUDING REMARKS 


Differential equations governing the thermal performance of charring abla- 
tors are derived. These equations are expanded into finite-difference form and 
have been programed for numerical solution on a digital computer. Numerical 
results compare favorably with available exact solutions. 

The required computer time presents a major problem in the numerical solu- 
tion of these equations. The governing equations have been examined and a num- 
ber of approximations which reduce the computer time are introduced. The con- 
ditions under which these simplifications should be used and the error involved 
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in their use are discussed. The computer program based on the equations pre- 
sented herein has been found to provide a practical basis for heat-shield 
design studies. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., April 1 , 1965* 
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APPENDIX A 


EFFECTS OF MASS INJECTION ON AERODYNAMIC HEATING 


The aerodynamic heat input consists of both convective heating and radia- 
tive heating. The material that is liberated by erosion of the outer surface 
and pyrolysis at the interface enters the boundary layer and causes a signifi- 
cant reduction in convective heating. It has been shown (refs. 20, 21, and 22) 
that, for moderate mass-transfer rates, the effectiveness of mass transfer in 
blocking convective heating is approximately a linear function of the enthalpy 
in the flow outside the boundary layer. It is further shown in references 22 
and 2 J, however, that this approximation may lead to serious error if the sur- 
face is exposed to a net radiant heat Input such as that which a typical vehi- 
cle would experience from reentry at parabolic or hyperbolic velocities. 

Blocking effectiveness as a function of a mass-transfer- rate parameter 
mh e /qQ ^where m = a c m c + Opip) is shown in figure 7* The exact solution was 

obtained from the boundary- layer solutions for air-to-air injection (ref. 21). 
It can be seen from figure 7 that the usual linear approximation for a laminar 
boundary layer is not valid at higher values of the mass-transfer-rate param- 
eter. The effect of radiant heating is to increase m and this increase 
results in operation at higher values of mh e /q£. The second-degree approxi- 
mation was developed by fitting a curve at values of mh^q^, of 0, 1.0, 
and 2 . 5 . 

Very little experimental data is available at high mass-transfer rates, 
and in conservative calculations it is desirable to specify a minimum value of 
^C,net 

greater than 0 . 


This approach is certainly 
in order if boundary- layer 
separation is found to 
occur at high mass- 
injection rates. In this 
analysis, a minimum value 
of the q- ratio (0.0k) was 
specified where the second- 
degree approximation departs 
from the theoretical curve. 

.2 

The equation for the 
second-degree approximation 
is as follows: 
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APPENDIX A 



The constants ac and. ctp are used to correct equation (Al) for the differ- 
ence between the molecular weight of the boundary- layer gas and the molecular 
weight of the injected gas. The constant a c must also be corrected for that 
part of the char that is removed mechanically, that is, for the char that is 
removed but not sublimed. A similar approach might be used to correct for 
turbulent flow. The effects of molecular weight and turbulent flow are dis- 
cussed in reference 2k. An equation of the form of equation (Al) is used to 
correlate experimental data in reference 18 by using the correction of refer- 
ence 2k for molecular weight effects. The equations have been formulated to 
provide the option of using either the linear function of enthalpy (ablation 
theory) or a second-degree approximation to the actual blocking effectiveness 
(transpiration theory). 
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APPENDIX B 


FINITE-DIFFERENCE EQUATIONS 


The system of partial differential equations to he solved is extremely 
complicated (nonlinear with variable coefficients) and solutions in analytic 
form can hardly be expected. The systems can be put into a form suitable for 
numerical calculation by deriving the equations in finite-difference form and 
then solving them by an iteration procedure on a high-speed digital computer. 

For this type of procedure it is convenient to space the finite-difference 
stations at equal intervals throughout a given layer, as shown in figure 2. 

The first (char) layer contains i stations including one at the front and one 
at the back surface of the first layer. The second layer contains j stations 
including one at the back surface of the second layer. The third layer (if 
desired) contains m stations including one at the back surface. 


The distances between stations are therefore 


l 



V = 



it 



(Bl) 


for the three layers whose thicknesses are x, x 1 , and x", respectively. In 
this appendix, the finite-difference form of the general equations (eqs. (24a), 
(24b), and (3)) are presented first and then the finite-difference form of the 
boundary equations (eqs. ( 25 a), ( 26 b), ( 27 ), and ( 29 )) are given. 


Interior Stations 


Char layer (2 ^ n ^ i - l).- The finite-difference form of equation (24a) 


i_ A/* d§\ 

x2B|\ <%) 


+ 


1 Se 


m c c p 


+ mpCp + | 



P 



+ F = 


„ ae 

pCp &t 


(B2) 


is obtained as follows: The second derivative (first term of eq. (B2)) is 

obtained by a Taylor series expansion at station n evaluated at the points 

^n - i) and ^n + ij. Thus, 
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Be\ = k /W\ _ A A 4 8e\ + (^) 2 &L4 Be\ _ (^) 3 £_4 
^1/ i \^i/n 2 d|\ B|/ n 8 S| 2 \ 8 | ) n 48 d|3\ 8|/ n 

n ' 2 


> (B3) 


Jte\ .. Jte) &t b_L Se^ + (^t) 2 £/ k 8e\ + (^) 3 8 3 L 3e\ 
\Slj 1 Wn 2 ■ %[ bil 8 ^2\ S| J n 48 S|5 \ ai/n 

o ^ 


These equations are then solved for the second derivative and yield 


i A4 ^ 


x 8|l L i 


ifk M 


*'»+§ 


Ak 


(Ag) 2 8^/ *e\ 

24x *Jn 


* n-i 

2 


(B4) 


or 


1 A( k |i 


T n+1 “ T n 


T n " T n-1 


X 


2 B| r B| L ~ l lfn,n+l - k n-l,n z 


O 

which is correct to order 2 * Note that x A| = 2. 

The finite-difference form of the first partial derivative in the second 
term of equation (B2) is obtained by a Taylor series expansion about station n 
evaluated at n - 1, n + 1, and n + 2. Thus 

e e - a(^) + (^.i 2 / sLi'i _ te) k (£e\ . 

n -l n *UA> ^VA, 6 yl * la ••• 


n ^[Zik 2 W 2 j n 


e n+l “ 9 


W - e„ ♦ *e(|) n - ^(0 


!&£(£*) + &£(£$ ) + . . 


« WL 2k W 


n 


> (®5) 


(2A| ) 5 /dV 


n 


V^ 3 /n 


(2ag)YdV\ . 

“Vi?y n 
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This set of equations is solved for the first derivative and yields 

* | = ir( 6T "l - 5T “ - <*> 

Therefore, the expression 

^■(6 T n +1 - 3T n - 2T n _l - T n+2 ) 

is an approximate derivative at station n (2 ^ n = i - 
order l ^> . 

At station n = i - 1, the series is evaluated at n - 2, n - 1, and 
n + 1 since there is a discontinuity at station i. This procedure gives 


(B7) 


2) correct to terms of 



= i(2Ti + JT,.! - 


6 T i-2 + T i-3) 


which is also correct to terms of order 


zl 


(b8) 


Substituting equations (B^) and (B7) into equation (B2) gives the finite- 
difference form of equation (B2) as 


•^n-1 ~ ^n 


■^n “ ^n+1 


k n-l,n i ~ k n,n+l j 
+ 7 (^n+1 “ 3T n - 2T n _q - T n+2 ^ 


. . - n - 1[ ™p p • ] 

m c c p + inpcp + - m c lc p 


AT. 

+ ZF = pepZ — 


n 


(2 ^ n < i - 2) (B9) 


and by using equation (b 8) instead of equation (B7)> the equation is 
T n-1 - T. 


k n-l, n 


L n 


- k n,n+l 


T n " T n+1 


+ 7 ( 2T n+l + 3T n - 6T n _q + T n _ 2 ^ 


n-1 V . ) 
m c c p + mpCp + I _ T l- r __ - m c l C p 


AT 

+ » - p= P J s 


n 


(n = i - 1) (BIO) 
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Uncharred material ( i + 1 ^ n ^ i + ,j - l) .- The transformed equation for 
the second layer (eq.. (24b) ) is 


1 d Lx c>0 *\ 1 ™P P ’ C P 


ae' 


jLf f (*’ w) + ^ ^ (1 ' °f~ = p ’^ tr 


(Bll) 


The second partial derivative is found as in the previous section and is 
given by the following finite-difference expression (correct to order 1%) • 


i _a/ k * M 

( x t )2 \ n 


T n-1 - T n 


T n “ T n+1 ] 


~~ ^ » rhi-1, n ^ t " -^n^ n+1 ^ i 


(B12) 


and as in the previous section, the first derivative is obtained by a Taylor 
series expansion about station n 


(~r = gp^n+l - 3T n - 2T n _! - T n+2 ) (i + 1 ^ n ^ i + J - 2) 


(B13) 


and 


(£ iLj.i ■ + 3T ° - * m 


(n = i + j - 1) (Bl4) 


These expressions are again correct to order 7 . 

This procedure gives the finite-difference form of equation (Bll) as 


, i T n-1 T n . i T n T n-1 
K n-1, n 2 ' - n +l 2 ’ 


+ |( 6T n+l " 3 T n - 2T n _! - T n+2 ) 


i + J 


n V ,C P = p > c >2* 

p' - p P V 1 At 


(i+l<|n<i + j-2) (B15) 
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and as 


T n-1 “ T n 


n-l,n z . 


- k ' 


T n “ T n+1 


n, n+1 i * 


+ |(2Tn+l + 3T n - 6T n _! + 1 n _ 2 ) - 1 p — = l 


iii ^n 

P P At 


(n = i + j - 1) (Bl6) 

Insulation (i + ,j + 1 ^ n ^ i + ,j + m- l) .- The differential equation for 
the third layer is (eq. ( 5 )) 


ST P °P St 


(B17) 


The second derivative is given by the expression (correct to order l^) 


JL/lr" 1 = V 

ay l )n " n ' 1 ’ n 


^n-1 “ ^n , " ^n “ ^n+1 
l " “ K n,n-1 z " 


(B18) 


Therefore, the difference equation is. 


^n-1 “ -^n ^n “ ^n+1 _ , ^-n 

K n-l,n Z " “ K n,n+1 z " 1 p °p ^t 


(B19) 


Boundary Stations 


First station (£ = 0) .- Some difficulty is encountered in finding an 
acceptable finite-difference equation for the first station. The first approach 
that is used to obtain such an equation is relatively straightforward. First, 
the second partial derivative of equation (B2) is expanded in a Taylor series 

Afc 

about the point £ = 0 and evaluated at | 


1 

X 2 



= 2 
l 


%2 




(B20) 
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Substituting equation (B20) into equation (B2) gives 


k gg_- . iL |i 

1,2 1 x v ^4=0 


^(f)l=o^ cCp + ^ 


+ F = pc p (E21) 


The transformed equation for the boundary condition at the first station 
(eq. ( 25 a)) is 


Wo - ae i° k - i if + *<= 


s ( e - T l)( H c + Ah 0 ) - ih c 


(B22) 


where q aero is given by. equation ( 25 b). 


The exact value of (— from equation (B22) is then substituted 

\ x oi/i=o 

into equation (B2l) to yield 


%2 


- UA ♦ Ac 


s(e - T^fHo + ih c ) - Ah c 


^aero/ 


1 - 4( A c c P + Vp) 


7 7 ^1 

+ 2 F " PC P 2 At 


(B23) 


The difficulty encountered in attempting to use this equation is due to the 
term 1 - ~^m c Cp + ipCpj . This term must be approximately equal to 1 since l 

is an approximation to an infinitesimal. However, since k is typically of 
the order 10 ”^ or smaller, l must be taken so small that the computer time 
required for an acceptable solution is unreasonable. 


To alleviate this problem, a different expression for the second first 
derivative of equation (B2l) is used. This expression was found by expanding 
the temperature at station 1 in a Taylor series evaluated at stations 2, 3 , 
and 4 which gives 


l(f) 0 ‘ 6t(' 11Ti + l8T2 ' 9T J + 


(B2U) 
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Hence the finite-difference equation for station 1 becomes 


1,2 “ + q aero " CTe l T i " “ c 


5(Ti - T^^Hc + Ah c ) - Ah c 


+ i(“ 11T l + l8T2 ' ^5 + ^(“cCp + ipc p ) + \ F = pc p | (B25) 

Note that the exact expression for ( — ^ ) from equation (B22) is still 

r ^ 4=0 


used for the first first derivative. 


Interface ( £ = 1, £ = 0) .- The transformed equation at the interface 

boundary ^ = 1, £ = 0 is (eq. ( 26 b)) 


k de\ 

, x 


= nip Ati p 


k_|_ del 

*’ k=o 


(B26) 


The finite-difference form of equation (B27) is found by a simultaneous 
solution of equations (B2) at £ = 1, (Bll) at £ = 0, and (B26). In this way, 
the interface equation is satisfied both at the interface and also immediately 
on either side of the interface. To do this, equation (B2l) is rewritten for 
£ = 1. Thus, 


"V X 




Ti 


- Tj.i 

l 


+ 


i /5e\ 

Sx\9| 



+ 



p c p * 


Z ^1 


2 At 


Then by using a Taylor’s series expansion about station i 
i - 2, and i - 3.? equation (B28) can be written as 


/k de\ _ t T i - T i-1 *p4 c P p V 

V* ¥j|.i " kl - 1 - 1 I + isl'P + p-r^J( 


y\rp < 

+ 9Ti_ 2 - 2Ti_ 3 ) + | F - p Cp | 



, evaluated at 
UTi - iSTi.i 


(B27) 

i - 1, 


(B28) 
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A similar procedure using equation (Bll) at £ = 0 gives 


k’ 8e 




sr h «> " ' ki ' i+1 ~ ~ 


T i+1 “ T i \>( C P P 


^=0 


12 \P - PI 


tV- -llTj. + l8T i+1 


\ i 7 1 AT-s 

- 9T 1+2 + 2T i+J ) + p'o p U _i 


(B29) 


Finally, substitution of equations (B28) and (B29) into equation (B26) gives 


T i+1 “ T i 


ki,i+l p *1-1,1 


^ - Tj..! , ^ ^ P C P 

+ “P K + 


p' - P 


i( 11T l - 


i8t 


i-l 


+ 9Ti_2 - 2Ti_ 3 ) 


c pP' 


p' - P 


i(-ll Tl t 18 t 1+1 - 9T 1+2 + 2T 1+3 ) 


- Ah T 


+ I F * ( pc P I + p,c P 


(B30) 


Back surface (2 layers, £ = 1).- For a two-layer system, the transformed 
boundary equation is (eq. ( 27 )) 

- ^(Ir ) ?=1 * lr + s ( e ' - 


+ TO i +J 


(e 1 )^ - 


(B31) 


The first term of equation (B3l) is determined by evaluating equation (Bll) 
at £ = 1; thus. 


k' de' 


= k' T i+j ~ T i+j-l , V ^i+j 

\x- 8£ J^ =1 i+j-l.i+J l' P P 2 At 

Substituting equation (B32) into equation (B3l) gives 

^1^10 . « 1+j (4 +J _ 4 


(B32) 


s(t 1+J - T 1+ j)»f “ (p S T + °i+j) 


At 


(B33) 
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Interface (3 layers, C = l) .- For a three-layer system, the transformed 
b oundary equation at £ = 1 is (e q . ( 28 b ) ) 


'k> 8e'\ _ p 8e ' 

F arJc-i 1+J a" 


L " z £\ 


(B5^) 


The first term of equation (B^k) is given by equation (Bj>2). The last 
term of equation (Bjk) is evaluated by expanding the second derivative of equa- 
tion (Bl8). Thus, 


- k 


89 " 

Sy 


x 0 +x 0 


- 

“ K i+j,i+j+l 


T i+j+l " T i+j 


+ p' 


2 


l" ^i+j 


At 


(B35) 


Substituting equations (BJ2) and (B35) into equation (B3*0 gives 


k i+j-i,i+j 


T i+j-l ~ T i+j 


T 


c i+j,i+j+l 


i+j " ^i+j+l 


, « ll 

i l r» [ 


u » l 


= \P ,c p 2 + P G p “o" + c i+j 


AT 


i+j 


At 


(B36) 


Back surface (3 layers, n - i + j + m) .- For a three-layer system, the 
back-surface boundary condition is (eq. ( 29 ) ) 


§Tr = -i+J+m §f- + s(e" - T 1+J+m )«f ihf + «i + j« 


(e”) 4 - Tg 


(B37) 


Equation (BJ7) is solved simultaneously with equation (B17), with the 
previously used approximation to the second derivative, to yield. 


a+j+m-1, i+j+m 


T i+j+m-l " T i+j+m 


- - 4) + s ( T i + J« - T i+J«.)«f + (pSt + C i+J+r.)^% !ia 


(BJ8) 
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SIMPLIFICATIONS TO REDUCE COMPUTING TIME 


Computing the Time Interval 


The primary difficulty encountered in using the program is the time 
required to obtain solutions on the digital computer. The equations are inte- 
grated stepwise over time, and the maximum time interval for which the solution 
is stable with the explicit formulation used is (ref. 25 ), 


At 


pepZ 2 

2 k 


where the properties of the layer giving the minimum At must be used. 

Theoretically, the initial value of Z may be zero in the physical sense; 
however, a finite initial char thickness is required for the calculating pro- 
cedure used in this program. The initial value of the char thickness which is 
usually chosen is extremely small. As the thickness of the char layer increases 
due to the pyrolysis at the interface, the preceding stability equation allows 
the digital computer to use a larger time interval in its stepwise calculations. 
Likewise, when the second material pyrolizes to such an extent that its Z 
becomes very small, the preceding equation now decreases the time interval so 
that the finite-difference equations remain stable. This time interval is only 
allowed to change whenever the digital computer can recognize this change and 
continue to print out answers at the specified print-out times. Also the com- 
puter recognizes specified minimum values for both x and x 1 so that com- 
puter time is not wasted when either x or x f approaches zero. 


Options for Thin-Layer Configurations 

In some cases, it is necessary to run calculations where the stability 
time interval At becomes very small. Therefore, to decrease computing time 
by at least one-half, the computer recognizes input quantities b and b ! . 
When x ^ b and/or x 1 ^ b ! , the interior stations are no longer considered 
and the entire layer becomes two half elements of thickness x/2 and/or x'/2- 
In this procedure the following finite-difference equations are utilized. 


When x ^ b, the equations at the first station (eq. (B25)) become 
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The equations for the interior stations in the char layer are no longer con- 
sidered when x ^ b. The equation at the interface (eqs. (B30)) is as follows: 

When x ^ b and x * > b ' 


2k 


i,i+l 


T i+1 " 

V 


2k- 




T i ~ T 1 
x 





i , 

Cpp’ 


|(-ll T i + l8T i+l - 9T i+2 + 2T i+5 )l - 2AhJl = (pc p x + P’c^Z')^ (C2) 


P* - p| 

When x ^ b and x ' ^ b ' 


2k. 


i T i+j “ T i 
i,i+j x' 


Ti - T-l 


2^1, i ' x " + xF + m p 


C— + 


pc T 


P ' p* - p '\ Ti " Tl 
= (pcpx + P' c p x ')Sr 


c p p ' / \ 

+ - Ti ) - ^ 


When x > b and x f < b 


) 


(C3) 


2k 


T i+j " T i T i - T i-1 

i,i+j x . “ ^i-lji i 


pc. 


+ 2F + Ms + ^ K - l8T i-i 


p' - p. 


+ ^1-2 - ^i-s)] + ^r7( T i+j - T 0 - - (pV + cS*')® 1 

When x* ^ b T , the interior stations in the second layer are not con- 
sidered. The equations at the back of the second layer (eq. (B33)) for two 
layers, (eq. (B 38 ) for three layers) are as follows: 


(Ok) 


When x < b 1 , equation (B33) becomes 


k i,i+J 


T i - 


T 


i+j 


x 1 


ae 


i+j 





Ah f 




(C5) 
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and equation (B38) becomes 


, , Ti - T i+J „ T i+j - T i+j+1 

k i,i+j x » " K i+j,i+j+l i” 




P c p 




AT 


i+d 


At 


(c 6) 

The computer uses the standard equations whenever x > b, and/or x’ > b*. 
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